Introduction

To organise input data for segmented regression of age-standardised mortality rates. * Load packages * Load data * Ensure correctly labelled * Convert dates to identifiable quarters * Visualise trends

Load packages

#install.packages("pacman")
pacman::p_load(
  tidyverse,
  segmented,
  plotly
               )

Load data and assign to object “data”

data <- read_csv("Data/QuarterlyASMR1990_2018.csv")
Parsed with column specification:
cols(
  period = col_double(),
  sex = col_character(),
  allageesprate = col_double(),
  allagestderr = col_double(),
  allagelcli = col_double(),
  allageucli = col_double(),
  under75esprate = col_double(),
  under75stderr = col_double(),
  under75lcli = col_double(),
  under75ucli = col_double(),
  totdeaths = col_double()
)

Look at data contents

data
glimpse(data)
Observations: 342
Variables: 11
$ period         <dbl> 19901, 19901, 19901, 19902, 19902, 19902, 19903, 19903, 19903, 19904, 19904, 19904, 19911, 19911, 19911, 1...
$ sex            <chr> "F", "M", "P", "F", "M", "P", "F", "M", "P", "F", "M", "P", "F", "M", "P", "F", "M", "P", "F", "M", "P", "...
$ allageesprate  <dbl> 1500.400, 2287.741, 1796.271, 1499.669, 2271.778, 1790.583, 1485.507, 2252.730, 1774.039, 1386.232, 2114.4...
$ allagestderr   <dbl> 7.788046, 13.950236, 6.854256, 7.774313, 13.869048, 6.830757, 7.730787, 13.786253, 6.788831, 7.495202, 13....
$ allagelcli     <dbl> 1485.136, 2260.399, 1782.837, 1484.431, 2244.595, 1777.194, 1470.355, 2225.709, 1760.733, 1371.541, 2088.2...
$ allageucli     <dbl> 1515.665, 2315.084, 1809.706, 1514.907, 2298.962, 1803.971, 1500.660, 2279.751, 1787.345, 1400.923, 2140.6...
$ under75esprate <dbl> 589.2652, 1010.2307, 777.9763, 590.4938, 1002.3490, 775.1230, 588.0220, 993.3030, 770.2041, 563.1330, 956....
$ under75stderr  <dbl> 4.994766, 7.222368, 4.242511, 4.998890, 7.183408, 4.230873, 4.985508, 7.136414, 4.212956, 4.878688, 6.9951...
$ under75lcli    <dbl> 579.4755, 996.0749, 769.6609, 580.6960, 988.2696, 766.8305, 578.2504, 979.3157, 761.9467, 553.5707, 942.94...
$ under75ucli    <dbl> 599.0550, 1024.3866, 786.2916, 600.2917, 1016.4285, 783.4155, 597.7936, 1007.2904, 778.4615, 572.6952, 970...
$ totdeaths      <dbl> 34163, 31360, 65523, 34216, 31247, 65463, 33997, 31069, 65066, 31910, 29617, 61527, 31664, 29208, 60872, 3...

Separate combined dates to years and quarters

data %>%
  separate(period, into = c("year","quarter"), sep = 4) %>% 
  mutate(
    year = as.double(year),
    quarter = as.double(quarter)
  ) %>% 
  mutate(
    time = year + (quarter / 4) - 0.125
  ) %>% 
  select(time, sex, totdeaths) -> tidied_totdeaths

Visualise number of deaths by sex

tidied_totdeaths %>%
  ggplot(
    aes(
      x = time, 
      y = totdeaths,
      color = sex
    )
  ) +
  geom_line() +
  coord_cartesian(ylim = c(0, 70000)) -> alldeathsplot

Visualise with plotly

alldeathsplot %>%
  ggplotly()

Visualise deaths men and women only

tidied_totdeaths %>%
  filter(
    sex != "P"
  ) %>%
  ggplot(
    aes(
      x = time, 
      y = totdeaths,
      color = sex
    )
  ) +
  geom_line() 

Tidy data - years and quarters into doubles; keep only all age, assign to object “tidied_allagerates”

data %>%
  separate(period, into = c("year","quarter"), sep = 4) %>% 
  mutate(
    year = as.double(year),
    quarter = as.double(quarter)
  ) %>% 
  mutate(
    time = year + (quarter / 4) - 0.125
  ) %>% 
  select(time, sex, allageesprate, allagelcli, allageucli) -> tidied_allagerates 

Visualise ASMR with ribbon

tidied_allagerates %>%
  filter(
    sex != "P"
  ) %>%
  ggplot(
    aes(
      x = time, 
      group = sex
    )
  ) +
  geom_ribbon(
    aes(
      ymin = allagelcli,
      ymax = allageucli,
      fill = sex
      )
  ) + 
    geom_line(
    aes(y = allageesprate)
  ) + 
  coord_cartesian(
    ylim = c(0, 3000)
  )  

Segmented regression

Generate linear model for both sexes

tidied_allagerates %>% 
  rename(rate = allageesprate) %>%
  filter(sex == "M") -> allagemale
  lm(rate ~ time, data = allagemale) -> mod_male
tidied_allagerates %>%
  rename(rate = allageesprate) %>%
  filter(sex == "F") -> allagefemale
  lm(rate ~ time, data = allagefemale) -> mod_female

View linear model

summary (mod_male)

Call:
lm(formula = rate ~ time, data = allagemale)

Residuals:
   Min     1Q Median     3Q    Max 
-68.73 -38.87 -16.80  36.84 142.55 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 69951.244   1164.509   60.07   <2e-16 ***
time          -34.056      0.581  -58.62   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 51.04 on 112 degrees of freedom
Multiple R-squared:  0.9684,    Adjusted R-squared:  0.9681 
F-statistic:  3436 on 1 and 112 DF,  p-value: < 2.2e-16
summary (mod_female)

Call:
lm(formula = rate ~ time, data = allagefemale)

Residuals:
    Min      1Q  Median      3Q     Max 
-56.519 -24.836  -3.137  17.272  89.669 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 34630.7949   734.8493   47.13   <2e-16 ***
time          -16.6907     0.3666  -45.52   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 32.21 on 112 degrees of freedom
Multiple R-squared:  0.9487,    Adjusted R-squared:  0.9483 
F-statistic:  2072 on 1 and 112 DF,  p-value: < 2.2e-16

Run one break segmented regression for both sexes and view

segmented(mod_male, seg.Z = ~time, psi = 1994.125) -> maleonebreak
segmented(mod_female, seg.Z = ~time, psi = 1994.125) -> femaleonebreak
summary(maleonebreak)

    ***Regression Model with Segmented Relationship(s)***

Call: 
segmented.lm(obj = mod_male, seg.Z = ~time, psi = 1994.125)

Estimated Break-Point(s):
               Est. St.Err
psi1.time 2013.281  0.474

Meaningful coefficients of the linear terms:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 76521.8711  1137.0750  67.297   <2e-16 ***
time          -37.3435     0.5681 -65.737   <2e-16 ***
U1.time        37.2230     5.3303   6.983       NA    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 36.77 on 110 degrees of freedom
Multiple R-Squared: 0.9839,  Adjusted R-squared: 0.9835 

Convergence attained in 2 iter. (rel. change 5.58e-07)
summary(femaleonebreak)

    ***Regression Model with Segmented Relationship(s)***

Call: 
segmented.lm(obj = mod_female, seg.Z = ~time, psi = 1994.125)

Estimated Break-Point(s):
               Est. St.Err
psi1.time 2014.375  0.576

Meaningful coefficients of the linear terms:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 37362.1905   793.9718  47.057   <2e-16 ***
time          -18.0571     0.3966 -45.534   <2e-16 ***
U1.time        24.0840     5.4286   4.437       NA    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 27.34 on 110 degrees of freedom
Multiple R-Squared: 0.9637,  Adjusted R-squared: 0.9627 

Convergence attained in 1 iter. (rel. change -3.5754e-08)

Generate predicted values based on model, and plot alongside observed values

tidied_allagerates %>%
  filter(sex == "M") %>% 
  mutate(
    male_onebreak_pred = predict(
      maleonebreak, 
      newdata = tidied_allagerates %>% 
        filter(sex == "M") %>% 
        select(time)
      )
    ) -> malepredicted 
tidied_allagerates %>%
  filter(sex == "F") %>% 
  mutate(
    female_onebreak_pred = predict(
      femaleonebreak, 
      newdata = tidied_allagerates %>% 
        filter(sex == "F") %>% 
        select(time)
      )
    ) -> femalepredicted 
malepredicted %>%
  ggplot(
    aes(x = time)
  ) + 
  geom_point(
    aes(y = allageesprate)
  ) + 
  geom_line(
    aes(y = male_onebreak_pred)
  )

Merge male and female predicted values and assign to object “obs_pred_both”

femalepredicted %>% 
  select(
    time, sex, observed = allageesprate,
    lower = allagelcli, upper = allageucli,
    predicted = female_onebreak_pred
  ) %>% 
  bind_rows(
    malepredicted %>% 
      select(
        time, sex, observed = allageesprate,
        lower = allagelcli, upper = allageucli,
        predicted = male_onebreak_pred
      ) 
  ) -> obs_pred_both

Plot male and female observed and predicted

obs_pred_both %>%
  ggplot(
    aes(x = time, color = sex)
  ) + 
  geom_point(
    aes(y = observed)
  ) + 
  geom_line(
    aes(y = predicted)
  )

Generate and plot residuals based on one break model, x intercepts manually entered

obs_pred_both %>% 
  mutate(residual=observed-predicted) %>% 
  ggplot(
    aes(x = time, shape = sex, color=sex)
  ) + 
  geom_point(
    aes(y = residual)
  )+ 
  ylim(c(-100,100))+
  geom_hline(yintercept = 0)+
  geom_vline(xintercept = 2013.496)+
  geom_vline(xintercept = 2014.373, linetype= "dashed") 

Generate and plot residuals, plotting sexes on separate charts, x intercepts drawn direct from model

obs_pred_both %>% 
  mutate(residual=observed-predicted) %>% 
  ggplot(
    aes(x = time)
  ) + 
  geom_point(
    aes(y = residual)
  ) +
  facet_wrap(~sex) +
  ylim(c(-100,100))+
  geom_hline(yintercept = 0)+
  geom_vline(xintercept = maleonebreak[["psi"]][[2]])+
    #geom_vline(xintercept = 2013.496)+
  geom_vline(xintercept = femaleonebreak[["psi"]][[2]], linetype= "dashed") 

    #geom_vline(xintercept = 2014.373, linetype= "dashed") 

Two break segmented model

segmented(mod_male, seg.Z = ~time, psi = c(1994, 2002)) -> maletwobreaks
segmented(mod_female, seg.Z = ~time, psi = c(1994, 2002)) -> femaletwobreaks
LS0tDQp0aXRsZTogIkRhdGEgbWFuYWdlbWVudCBmb3Igc2VnbWVudGVkIHJlZ3Jlc3Npb24iDQpvdXRwdXQ6DQogIGh0bWxfbm90ZWJvb2s6IGRlZmF1bHQNCiAgaHRtbF9kb2N1bWVudDoNCiAgICBkZl9wcmludDogcGFnZWQNCiAgd29yZF9kb2N1bWVudDogZGVmYXVsdA0KLS0tDQojIEludHJvZHVjdGlvbg0KDQpUbyBvcmdhbmlzZSBpbnB1dCBkYXRhIGZvciBzZWdtZW50ZWQgcmVncmVzc2lvbiBvZiBhZ2Utc3RhbmRhcmRpc2VkIG1vcnRhbGl0eSByYXRlcy4gDQoqIExvYWQgcGFja2FnZXMNCiogTG9hZCBkYXRhDQoqIEVuc3VyZSBjb3JyZWN0bHkgbGFiZWxsZWQNCiogQ29udmVydCBkYXRlcyB0byBpZGVudGlmaWFibGUgcXVhcnRlcnMNCiogVmlzdWFsaXNlIHRyZW5kcw0KDQpgYGB7ciBzZXR1cCwgaW5jbHVkZT1GQUxTRSwgZWNobz1GQUxTRX0NCnJlcXVpcmUoImtuaXRyIikNCmtuaXRyOjpvcHRzX2tuaXQkc2V0KHJvb3QuZGlyID0gcnByb2pyb290OjpmaW5kX3JzdHVkaW9fcm9vdF9maWxlKCkpDQpgYGANCg0KIyMgTG9hZCBwYWNrYWdlcw0KDQpgYGB7cn0NCiNpbnN0YWxsLnBhY2thZ2VzKCJwYWNtYW4iKQ0KcGFjbWFuOjpwX2xvYWQoDQogIHRpZHl2ZXJzZSwNCiAgc2VnbWVudGVkLA0KICBwbG90bHkNCiAgICAgICAgICAgICAgICkNCmBgYA0KDQojIyBMb2FkIGRhdGEgYW5kIGFzc2lnbiB0byBvYmplY3QgImRhdGEiDQoNCmBgYHtyfQ0KZGF0YSA8LSByZWFkX2NzdigiRGF0YS9RdWFydGVybHlBU01SMTk5MF8yMDE4LmNzdiIpDQoNCmBgYA0KDQojIyBMb29rIGF0IGRhdGEgY29udGVudHMNCg0KYGBge3J9DQpkYXRhDQpnbGltcHNlKGRhdGEpDQpgYGANCg0KIyMgU2VwYXJhdGUgY29tYmluZWQgZGF0ZXMgdG8geWVhcnMgYW5kIHF1YXJ0ZXJzDQoNCmBgYHtyfQ0KZGF0YSAlPiUNCiAgc2VwYXJhdGUocGVyaW9kLCBpbnRvID0gYygieWVhciIsInF1YXJ0ZXIiKSwgc2VwID0gNCkgJT4lIA0KICBtdXRhdGUoDQogICAgeWVhciA9IGFzLmRvdWJsZSh5ZWFyKSwNCiAgICBxdWFydGVyID0gYXMuZG91YmxlKHF1YXJ0ZXIpDQogICkgJT4lIA0KICBtdXRhdGUoDQogICAgdGltZSA9IHllYXIgKyAocXVhcnRlciAvIDQpIC0gMC4xMjUNCiAgKSAlPiUgDQogIHNlbGVjdCh0aW1lLCBzZXgsIHRvdGRlYXRocykgLT4gdGlkaWVkX3RvdGRlYXRocw0KYGBgDQoNCg0KIyMgVmlzdWFsaXNlIG51bWJlciBvZiBkZWF0aHMgYnkgc2V4DQoNCmBgYHtyfQ0KdGlkaWVkX3RvdGRlYXRocyAlPiUNCiAgZ2dwbG90KA0KICAgIGFlcygNCiAgICAgIHggPSB0aW1lLCANCiAgICAgIHkgPSB0b3RkZWF0aHMsDQogICAgICBjb2xvciA9IHNleA0KICAgICkNCiAgKSArDQogIGdlb21fbGluZSgpICsNCiAgY29vcmRfY2FydGVzaWFuKHlsaW0gPSBjKDAsIDcwMDAwKSkgLT4gYWxsZGVhdGhzcGxvdA0KDQpgYGANCg0KIyMgVmlzdWFsaXNlIHdpdGggcGxvdGx5DQoNCmBgYHtyfQ0KYWxsZGVhdGhzcGxvdCAlPiUNCiAgZ2dwbG90bHkoKQ0KDQpgYGANCiMjIFZpc3VhbGlzZSBkZWF0aHMgbWVuIGFuZCB3b21lbiBvbmx5DQoNCmBgYHtyfQ0KdGlkaWVkX3RvdGRlYXRocyAlPiUNCiAgZmlsdGVyKA0KICAgIHNleCAhPSAiUCINCiAgKSAlPiUNCiAgZ2dwbG90KA0KICAgIGFlcygNCiAgICAgIHggPSB0aW1lLCANCiAgICAgIHkgPSB0b3RkZWF0aHMsDQogICAgICBjb2xvciA9IHNleA0KICAgICkNCiAgKSArDQogIGdlb21fbGluZSgpIA0KDQpgYGANCiMjIFRpZHkgZGF0YSAtIHllYXJzIGFuZCBxdWFydGVycyBpbnRvIGRvdWJsZXM7IGtlZXAgb25seSBhbGwgYWdlLCBhc3NpZ24gdG8gb2JqZWN0ICJ0aWRpZWRfYWxsYWdlcmF0ZXMiDQoNCmBgYHtyfQ0KZGF0YSAlPiUNCiAgc2VwYXJhdGUocGVyaW9kLCBpbnRvID0gYygieWVhciIsInF1YXJ0ZXIiKSwgc2VwID0gNCkgJT4lIA0KICBtdXRhdGUoDQogICAgeWVhciA9IGFzLmRvdWJsZSh5ZWFyKSwNCiAgICBxdWFydGVyID0gYXMuZG91YmxlKHF1YXJ0ZXIpDQogICkgJT4lIA0KICBtdXRhdGUoDQogICAgdGltZSA9IHllYXIgKyAocXVhcnRlciAvIDQpIC0gMC4xMjUNCiAgKSAlPiUgDQogIHNlbGVjdCh0aW1lLCBzZXgsIGFsbGFnZWVzcHJhdGUsIGFsbGFnZWxjbGksIGFsbGFnZXVjbGkpIC0+IHRpZGllZF9hbGxhZ2VyYXRlcyANCmBgYA0KDQojI1Zpc3VhbGlzZSBBU01SIHdpdGggcmliYm9uDQoNCmBgYHtyfQ0KdGlkaWVkX2FsbGFnZXJhdGVzICU+JQ0KICBmaWx0ZXIoDQogICAgc2V4ICE9ICJQIg0KICApICU+JQ0KICBnZ3Bsb3QoDQogICAgYWVzKA0KICAgICAgeCA9IHRpbWUsIA0KICAgICAgZ3JvdXAgPSBzZXgNCiAgICApDQogICkgKw0KDQogIGdlb21fcmliYm9uKA0KICAgIGFlcygNCiAgICAgIHltaW4gPSBhbGxhZ2VsY2xpLA0KICAgICAgeW1heCA9IGFsbGFnZXVjbGksDQogICAgICBmaWxsID0gc2V4DQogICAgICApDQogICkgKyANCiAgICBnZW9tX2xpbmUoDQogICAgYWVzKHkgPSBhbGxhZ2Vlc3ByYXRlKQ0KICApICsgDQogIGNvb3JkX2NhcnRlc2lhbigNCiAgICB5bGltID0gYygwLCAzMDAwKQ0KICApICANCmBgYA0KDQojIFNlZ21lbnRlZCByZWdyZXNzaW9uIA0KIyNHZW5lcmF0ZSBsaW5lYXIgbW9kZWwgZm9yIGJvdGggc2V4ZXMNCg0KYGBge3J9DQp0aWRpZWRfYWxsYWdlcmF0ZXMgJT4lIA0KICByZW5hbWUocmF0ZSA9IGFsbGFnZWVzcHJhdGUpICU+JQ0KICBmaWx0ZXIoc2V4ID09ICJNIikgLT4gYWxsYWdlbWFsZQ0KICBsbShyYXRlIH4gdGltZSwgZGF0YSA9IGFsbGFnZW1hbGUpIC0+IG1vZF9tYWxlDQoNCnRpZGllZF9hbGxhZ2VyYXRlcyAlPiUNCiAgcmVuYW1lKHJhdGUgPSBhbGxhZ2Vlc3ByYXRlKSAlPiUNCiAgZmlsdGVyKHNleCA9PSAiRiIpIC0+IGFsbGFnZWZlbWFsZQ0KICBsbShyYXRlIH4gdGltZSwgZGF0YSA9IGFsbGFnZWZlbWFsZSkgLT4gbW9kX2ZlbWFsZQ0KDQpgYGANCg0KIyNWaWV3IGxpbmVhciBtb2RlbCANCg0KYGBge3J9DQpzdW1tYXJ5IChtb2RfbWFsZSkNCnN1bW1hcnkgKG1vZF9mZW1hbGUpDQoNCmBgYA0KDQojI1J1biBvbmUgYnJlYWsgc2VnbWVudGVkIHJlZ3Jlc3Npb24gZm9yIGJvdGggc2V4ZXMgYW5kIHZpZXcNCg0KYGBge3J9DQpzZWdtZW50ZWQobW9kX21hbGUsIHNlZy5aID0gfnRpbWUsIHBzaSA9IDE5OTQuMTI1KSAtPiBtYWxlb25lYnJlYWsNCnNlZ21lbnRlZChtb2RfZmVtYWxlLCBzZWcuWiA9IH50aW1lLCBwc2kgPSAxOTk0LjEyNSkgLT4gZmVtYWxlb25lYnJlYWsNCnN1bW1hcnkobWFsZW9uZWJyZWFrKQ0Kc3VtbWFyeShmZW1hbGVvbmVicmVhaykNCg0KYGBgDQoNCiMjR2VuZXJhdGUgcHJlZGljdGVkIHZhbHVlcyBiYXNlZCBvbiBtb2RlbCwgYW5kIHBsb3QgYWxvbmdzaWRlIG9ic2VydmVkIHZhbHVlcw0KDQpgYGB7cn0NCnRpZGllZF9hbGxhZ2VyYXRlcyAlPiUNCiAgZmlsdGVyKHNleCA9PSAiTSIpICU+JSANCiAgbXV0YXRlKA0KICAgIG1hbGVfb25lYnJlYWtfcHJlZCA9IHByZWRpY3QoDQogICAgICBtYWxlb25lYnJlYWssIA0KICAgICAgbmV3ZGF0YSA9IHRpZGllZF9hbGxhZ2VyYXRlcyAlPiUgDQogICAgICAgIGZpbHRlcihzZXggPT0gIk0iKSAlPiUgDQogICAgICAgIHNlbGVjdCh0aW1lKQ0KICAgICAgKQ0KICAgICkgLT4gbWFsZXByZWRpY3RlZCANCg0KdGlkaWVkX2FsbGFnZXJhdGVzICU+JQ0KICBmaWx0ZXIoc2V4ID09ICJGIikgJT4lIA0KICBtdXRhdGUoDQogICAgZmVtYWxlX29uZWJyZWFrX3ByZWQgPSBwcmVkaWN0KA0KICAgICAgZmVtYWxlb25lYnJlYWssIA0KICAgICAgbmV3ZGF0YSA9IHRpZGllZF9hbGxhZ2VyYXRlcyAlPiUgDQogICAgICAgIGZpbHRlcihzZXggPT0gIkYiKSAlPiUgDQogICAgICAgIHNlbGVjdCh0aW1lKQ0KICAgICAgKQ0KICAgICkgLT4gZmVtYWxlcHJlZGljdGVkIA0KDQoNCm1hbGVwcmVkaWN0ZWQgJT4lDQogIGdncGxvdCgNCiAgICBhZXMoeCA9IHRpbWUpDQogICkgKyANCiAgZ2VvbV9wb2ludCgNCiAgICBhZXMoeSA9IGFsbGFnZWVzcHJhdGUpDQogICkgKyANCiAgZ2VvbV9saW5lKA0KICAgIGFlcyh5ID0gbWFsZV9vbmVicmVha19wcmVkKQ0KICApDQpgYGANCg0KIyNNZXJnZSBtYWxlIGFuZCBmZW1hbGUgcHJlZGljdGVkIHZhbHVlcyBhbmQgYXNzaWduIHRvIG9iamVjdCAib2JzX3ByZWRfYm90aCINCg0KYGBge3J9DQpmZW1hbGVwcmVkaWN0ZWQgJT4lIA0KICBzZWxlY3QoDQogICAgdGltZSwgc2V4LCBvYnNlcnZlZCA9IGFsbGFnZWVzcHJhdGUsDQogICAgbG93ZXIgPSBhbGxhZ2VsY2xpLCB1cHBlciA9IGFsbGFnZXVjbGksDQogICAgcHJlZGljdGVkID0gZmVtYWxlX29uZWJyZWFrX3ByZWQNCiAgKSAlPiUgDQogIGJpbmRfcm93cygNCiAgICBtYWxlcHJlZGljdGVkICU+JSANCiAgICAgIHNlbGVjdCgNCiAgICAgICAgdGltZSwgc2V4LCBvYnNlcnZlZCA9IGFsbGFnZWVzcHJhdGUsDQogICAgICAgIGxvd2VyID0gYWxsYWdlbGNsaSwgdXBwZXIgPSBhbGxhZ2V1Y2xpLA0KICAgICAgICBwcmVkaWN0ZWQgPSBtYWxlX29uZWJyZWFrX3ByZWQNCiAgICAgICkgDQogICkgLT4gb2JzX3ByZWRfYm90aA0KDQoNCmBgYA0KDQojI1Bsb3QgbWFsZSBhbmQgZmVtYWxlIG9ic2VydmVkIGFuZCBwcmVkaWN0ZWQNCg0KYGBge3J9DQpvYnNfcHJlZF9ib3RoICU+JQ0KICBnZ3Bsb3QoDQogICAgYWVzKHggPSB0aW1lLCBjb2xvciA9IHNleCkNCiAgKSArIA0KICBnZW9tX3BvaW50KA0KICAgIGFlcyh5ID0gb2JzZXJ2ZWQpDQogICkgKyANCiAgZ2VvbV9saW5lKA0KICAgIGFlcyh5ID0gcHJlZGljdGVkKQ0KICApDQpgYGANCg0KIyNHZW5lcmF0ZSBhbmQgcGxvdCByZXNpZHVhbHMgYmFzZWQgb24gb25lIGJyZWFrIG1vZGVsLCB4IGludGVyY2VwdHMgbWFudWFsbHkgZW50ZXJlZA0KDQpgYGB7cn0NCm9ic19wcmVkX2JvdGggJT4lIA0KICBtdXRhdGUocmVzaWR1YWw9b2JzZXJ2ZWQtcHJlZGljdGVkKSAlPiUgDQogIGdncGxvdCgNCiAgICBhZXMoeCA9IHRpbWUsIHNoYXBlID0gc2V4LCBjb2xvcj1zZXgpDQogICkgKyANCiAgZ2VvbV9wb2ludCgNCiAgICBhZXMoeSA9IHJlc2lkdWFsKQ0KICApKyANCiAgeWxpbShjKC0xMDAsMTAwKSkrDQogIGdlb21faGxpbmUoeWludGVyY2VwdCA9IDApKw0KICBnZW9tX3ZsaW5lKHhpbnRlcmNlcHQgPSAyMDEzLjQ5NikrDQogIGdlb21fdmxpbmUoeGludGVyY2VwdCA9IDIwMTQuMzczLCBsaW5ldHlwZT0gImRhc2hlZCIpIA0KDQpgYGANCg0KIyNHZW5lcmF0ZSBhbmQgcGxvdCByZXNpZHVhbHMsIHBsb3R0aW5nIHNleGVzIG9uIHNlcGFyYXRlIGNoYXJ0cywgeCBpbnRlcmNlcHRzIGRyYXduIGRpcmVjdCBmcm9tIG1vZGVsDQoNCmBgYHtyfQ0Kb2JzX3ByZWRfYm90aCAlPiUgDQogIG11dGF0ZShyZXNpZHVhbD1vYnNlcnZlZC1wcmVkaWN0ZWQpICU+JSANCiAgZ2dwbG90KA0KICAgIGFlcyh4ID0gdGltZSkNCiAgKSArIA0KICBnZW9tX3BvaW50KA0KICAgIGFlcyh5ID0gcmVzaWR1YWwpDQogICkgKw0KICBmYWNldF93cmFwKH5zZXgpICsNCiAgeWxpbShjKC0xMDAsMTAwKSkrDQogIGdlb21faGxpbmUoeWludGVyY2VwdCA9IDApKw0KICBnZW9tX3ZsaW5lKHhpbnRlcmNlcHQgPSBtYWxlb25lYnJlYWtbWyJwc2kiXV1bWzJdXSkrDQogICAgI2dlb21fdmxpbmUoeGludGVyY2VwdCA9IDIwMTMuNDk2KSsNCiAgZ2VvbV92bGluZSh4aW50ZXJjZXB0ID0gZmVtYWxlb25lYnJlYWtbWyJwc2kiXV1bWzJdXSwgbGluZXR5cGU9ICJkYXNoZWQiKSANCiAgICAjZ2VvbV92bGluZSh4aW50ZXJjZXB0ID0gMjAxNC4zNzMsIGxpbmV0eXBlPSAiZGFzaGVkIikgDQoNCg0KYGBgDQoNCiMjVHdvIGJyZWFrIHNlZ21lbnRlZCBtb2RlbCAgICANCg0KYGBge3J9DQpzZWdtZW50ZWQobW9kX21hbGUsIHNlZy5aID0gfnRpbWUsIHBzaSA9IGMoMTk5NCwgMjAwMikpIC0+IG1hbGV0d29icmVha3MNCnNlZ21lbnRlZChtb2RfZmVtYWxlLCBzZWcuWiA9IH50aW1lLCBwc2kgPSBjKDE5OTQsIDIwMDIpKSAtPiBmZW1hbGV0d29icmVha3MNCg0KYGBg